// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //


#include "potential.h"

std::string background::constant::record_atom_index() const { return ""; }
std::string background::constant::record_key() const { return "background const"; }
std::string background::constant::record_params() const 
{ std::ostringstream os; os << ios_param << p_u0; return os.str(); }
bool background::constant::read(const std::string & input,const group_t & group,const topology_t& topology,potential_list_t & potential_list)
{
	param_type p_u0;
	if(io::match("background const %1",input,&p_u0))
	{
		potential_list.insert(new constant(p_u0));
		return true;
	}
	return false;
}
bool background::constant_fit::read(const std::string & input,const group_t & group,const topology_t& topology,potential_fit_list_t & potential_list)
{
	io::prefix_mark_expression<real_type,false> p_u0;
	vector_index_type i_u0 = -1;
	std::vector<std::string> ie;
	if(io::match("background const %1",input,&p_u0))
	{
		if(~p_u0) { i_u0 = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_u0); }

		potential_list.insert(new constant_fit(*p_u0,i_u0));
		return true;
	}
	return false;
}
	


std::string bond::potential::record_atom_index() const 
{ std::ostringstream os; os << function::ait2ai(i) << "," << function::ait2ai(j); return os.str(); }

std::string bond::harm::record_key() const { return "bond harm"; }
std::string bond::harm::record_params() const 
{ std::ostringstream os; os << ios_param << p_k << ios_param << p_r0; return os.str(); }
bool bond::harm::read(const std::string & input,const group_t & group,const topology_t& topology,potential_list_t & potential_list)
{
	param_type p_k,p_r0;
	std::vector<std::string> ie;
	if(io::match("bond harm %1 %2 : %3",input,&p_k,&p_r0,&ie))
	{
		atom_index i,j;
		std::string si,sj;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2",',',*item,&i,&j)) 
		{	
			if(i!=j) potential_list.insert(new harm(p_k,p_r0,function::ai2ait(i),function::ai2ait(j)));
		} 
		else if (io::match("%1,%2",',',*item,&si,&sj))
		{
			const std::set<atom_index> & i_set = group[si]; //! \todo What if it does not exist.
			const std::set<atom_index> & j_set = group[sj];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			{
				potential_list.insert(new harm(p_k,p_r0,function::ai2ait(*it),function::ai2ait(*jt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}
bool bond::harm_fit::read(const std::string & input,const group_t & group,const topology_t& topology,potential_fit_list_t & potential_list)
{
	io::prefix_mark_expression<real_type,false> p_k,p_r0;
	vector_index_type i_k = -1,i_r0 = -1;
	std::vector<std::string> ie;
	if(io::match("bond harm %1 %2 : %3",input,&p_k,&p_r0,&ie))
	{
		if(~p_k) { i_k = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_k); }
		if(~p_r0) { i_r0 = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_r0); }

		atom_index i,j;
		std::string si,sj;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2",',',*item,&i,&j)) 
		{	
			if(i!=j) potential_list.insert(new harm_fit(*p_k,i_k,*p_r0,i_r0,function::ai2ait(i),function::ai2ait(j)));
		} 
		else if (io::match("%1,%2",',',*item,&si,&sj))
		{
			const std::set<atom_index> & i_set = group[si]; //! \todo What if it does not exist.
			const std::set<atom_index> & j_set = group[sj];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			{
				potential_list.insert(new harm_fit(*p_k,i_k,*p_r0,i_r0,function::ai2ait(*it),function::ai2ait(*jt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}
			


std::string bond::mors::record_key() const { return "bond mors"; }
std::string bond::mors::record_params() const 
{ std::ostringstream os; os << ios_param << p_E0 << ios_param << p_k << ios_param << p_r0; return os.str(); }
bool bond::mors::read(const std::string & input,const group_t & group,const topology_t& topology,potential_list_t & potential_list)
{
	param_type p_E0,p_k,p_r0;
	std::vector<std::string> ie;
	if(io::match("bond mors %1 %2 %3 : %4",input,&p_E0,&p_k,&p_r0,&ie))
	{
		atom_index i,j;
		std::string si,sj;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2",',',*item,&i,&j)) 
		{	
			if(i!=j) potential_list.insert(new mors(p_E0,p_k,p_r0,function::ai2ait(i),function::ai2ait(j)));
		} 
		else if (io::match("%1,%2",',',*item,&si,&sj))
		{
			const std::set<atom_index> & i_set = group[si]; //! \todo What if it does not exist.
			const std::set<atom_index> & j_set = group[sj];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			{
				potential_list.insert(new mors(p_E0,p_k,p_r0,function::ai2ait(*it),function::ai2ait(*jt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}
bool bond::mors_fit::read(const std::string & input,const group_t & group,const topology_t& topology,potential_fit_list_t & potential_list)
{
	io::prefix_mark_expression<real_type,false> p_E0,p_k,p_r0;
	vector_index_type i_E0 = -1,i_k = -1,i_r0 = -1;
	std::vector<std::string> ie;
	if(io::match("bond mors %1 %2 %3 : %4",input,&p_E0,&p_k,&p_r0,&ie))
	{
		if(~p_E0) { i_E0 = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_E0); }
		if(~p_k)  {  i_k = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_k);  }
		if(~p_r0) { i_r0 = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_r0); }


		atom_index i,j;
		std::string si,sj;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2",',',*item,&i,&j)) 
		{	
			if(i!=j) potential_list.insert(new mors_fit(*p_E0,i_E0,*p_k,i_k,*p_r0,i_r0,function::ai2ait(i),function::ai2ait(j)));
		} 
		else if (io::match("%1,%2",',',*item,&si,&sj))
		{
			const std::set<atom_index> & i_set = group[si]; //! \todo What if it does not exist.
			const std::set<atom_index> & j_set = group[sj];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			{
				potential_list.insert(new mors_fit(*p_E0,i_E0,*p_k,i_k,*p_r0,i_r0,function::ai2ait(*it),function::ai2ait(*jt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}
			




std::string angle::potential::record_atom_index() const 
{ std::ostringstream os; os << function::ait2ai(i) << "," << function::ait2ai(j)<< "," << function::ait2ai(k); return os.str(); }
		
std::string angle::harm::record_key() const { return "angle harm"; }
std::string angle::harm::record_params() const 
{ std::ostringstream os; os << ios_param << p_k << ios_param << function::rad2deg(p_phi0); return os.str(); }			
bool angle::harm::read(const std::string & input,const group_t & group,const topology_t& topology,potential_list_t & potential_list)
{
	param_type p_k,p_phi0;
	std::vector<std::string> ie;
	if(io::match("angle harm %1 %2 : %3",input,&p_k,&p_phi0,&ie))
	{
		p_phi0 = function::deg2rad(p_phi0);
		atom_index i,j,k;
		std::string si,sj,sk;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2,%3",',',*item,&i,&j,&k)) 
		{	
			if(i!=j) potential_list.insert(new harm(p_k,p_phi0,function::ai2ait(i),function::ai2ait(j),function::ai2ait(k)));
		} 
		else if (io::match("%1,%2,%3",',',*item,&si,&sj,&sk))
		{
			const std::set<atom_index> & i_set = group[si];
			const std::set<atom_index> & j_set = group[sj];
			const std::set<atom_index> & k_set = group[sk];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			for(std::set<atom_index>::const_iterator kt = k_set.begin(); kt != k_set.end(); kt++)
			if(*kt != *it)
			if(*kt != *jt)
			if(topology.path(*jt,*kt))
			{
				potential_list.insert(new harm(p_k,p_phi0,function::ai2ait(*it),function::ai2ait(*jt),function::ai2ait(*kt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}

bool angle::harm_fit::read(const std::string & input,const group_t & group,const topology_t& topology,potential_fit_list_t & potential_list )
{
	io::prefix_mark_expression<real_type,false> p_k,p_phi0;
	vector_index_type i_k = -1,i_phi0 = -1;
	std::vector<std::string> ie;
	if(io::match("angle harm %1 %2 : %3",input,&p_k,&p_phi0,&ie))
	{
		*p_phi0 = function::deg2rad(*p_phi0);

		if(~p_k) { i_k = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_k); }
		if(~p_phi0) { i_phi0 = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_phi0); }

		atom_index i,j,k;
		std::string si,sj,sk;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2,%3",',',*item,&i,&j,&k)) 
		{	
			if(i!=j) potential_list.insert(new harm_fit(*p_k,i_k,*p_phi0,i_phi0,function::ai2ait(i),function::ai2ait(j),function::ai2ait(k)));
		} 
		else if (io::match("%1,%2,%3",',',*item,&si,&sj,&sk))
		{
			const std::set<atom_index> & i_set = group[si];
			const std::set<atom_index> & j_set = group[sj];
			const std::set<atom_index> & k_set = group[sk];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			for(std::set<atom_index>::const_iterator kt = k_set.begin(); kt != k_set.end(); kt++)
			if(*kt != *it)
			if(*kt != *jt)
			if(topology.path(*jt,*kt))
			{
				potential_list.insert(new harm_fit(*p_k,i_k,*p_phi0,i_phi0,function::ai2ait(*it),function::ai2ait(*jt),function::ai2ait(*kt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}




std::string angle::hcos::record_key() const { return "angle hcos"; }
std::string angle::hcos::record_params() const 
{ std::ostringstream os; os << ios_param << p_k << ios_param << function::rad2deg(acos(p_cosphi0)); return os.str(); }			
bool angle::hcos::read(const std::string & input,const group_t & group,const topology_t& topology,potential_list_t & potential_list)
{
	param_type p_k,p_cosphi0;
	std::vector<std::string> ie;
	if(io::match("angle hcos %1 %2 : %3",input,&p_k,&p_cosphi0,&ie))  //! Parameter cosphi0 is given by phi0 in deg!
	{
		p_cosphi0 = std::cos(function::deg2rad(p_cosphi0));
		atom_index i,j,k;
		std::string si,sj,sk;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2,%3",',',*item,&i,&j,&k)) 
		{	
			if(i!=j) potential_list.insert(new hcos(p_k,p_cosphi0,function::ai2ait(i),function::ai2ait(j),function::ai2ait(k)));
		} 
		else if (io::match("%1,%2,%3",',',*item,&si,&sj,&sk))
		{
			const std::set<atom_index> & i_set = group[si];
			const std::set<atom_index> & j_set = group[sj];
			const std::set<atom_index> & k_set = group[sk];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			for(std::set<atom_index>::const_iterator kt = k_set.begin(); kt != k_set.end(); kt++)
			if(*kt != *it)
			if(*kt != *jt)
			if(topology.path(*jt,*kt))
			{
				potential_list.insert(new hcos(p_k,p_cosphi0,function::ai2ait(*it),function::ai2ait(*jt),function::ai2ait(*kt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}

bool angle::hcos_fit::read(const std::string & input,const group_t & group,const topology_t& topology,potential_fit_list_t & potential_list)
{
	io::prefix_mark_expression<real_type,false> p_k,p_cosphi0;
	vector_index_type i_k = -1,i_cosphi0 = -1;
	std::vector<std::string> ie;
	if(io::match("angle hcos %1 %2 : %3",input,&p_k,&p_cosphi0,&ie))
	{
		*p_cosphi0 = std::cos(function::deg2rad(*p_cosphi0));
                
		if(~p_k) { i_k = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_k); }
		if(~p_cosphi0) { i_cosphi0 = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_cosphi0); }

		atom_index i,j,k;
		std::string si,sj,sk;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2,%3",',',*item,&i,&j,&k)) 
		{	
			//! \Todo i!=j j!=k k!=i  for all angle cases.
			if(i!=j) potential_list.insert(new hcos_fit(*p_k,i_k,*p_cosphi0,i_cosphi0,function::ai2ait(i),function::ai2ait(j),function::ai2ait(k)));
		} 
		else if (io::match("%1,%2,%3",',',*item,&si,&sj,&sk))
		{
			const std::set<atom_index> & i_set = group[si];
			const std::set<atom_index> & j_set = group[sj];
			const std::set<atom_index> & k_set = group[sk];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			for(std::set<atom_index>::const_iterator kt = k_set.begin(); kt != k_set.end(); kt++)
			if(*kt != *it)
			if(*kt != *jt)
			if(topology.path(*jt,*kt))
			{
				potential_list.insert(new hcos_fit(*p_k,i_k,*p_cosphi0,i_cosphi0,function::ai2ait(*it),function::ai2ait(*jt),function::ai2ait(*kt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}



//!



std::string dihedral::potential::record_atom_index() const 
{ std::ostringstream os; os << function::ait2ai(i) << "," << function::ait2ai(j)<< "," << function::ait2ai(k) << "," << function::ait2ai(l); return os.str(); }
		
std::string dihedral::harm::record_key() const { return "dihedral harm"; }
std::string dihedral::harm::record_params() const 
{ std::ostringstream os; os << ios_param << p_k << ios_param << function::rad2deg(p_phi0); return os.str(); }			
bool dihedral::harm::read(const std::string & input,const group_t & group,const topology_t& topology,potential_list_t & potential_list)
{
	param_type p_k,p_phi0;
	std::vector<std::string> ie;
	if(io::match("dihedral harm %1 %2 : %3",input,&p_k,&p_phi0,&ie))
	{
		p_phi0 = function::deg2rad(p_phi0);
		atom_index i,j,k,l;
		std::string si,sj,sk,sl;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2,%3,%4",',',*item,&i,&j,&k,&l)) 
		{	
			if((i!=j) && (i!=k) && (i!=l)  &&   (j!=k) && (j!=l)   &&  (k!=l) ) 
			potential_list.insert(new harm(p_k,p_phi0,function::ai2ait(i),function::ai2ait(j),function::ai2ait(k),function::ai2ait(l)));
		} 
		else if (io::match("%1,%2,%3,%4",',',*item,&si,&sj,&sk,&sl))
		{
			const std::set<atom_index> & i_set = group[si];
			const std::set<atom_index> & j_set = group[sj];
			const std::set<atom_index> & k_set = group[sk];
			const std::set<atom_index> & l_set = group[sl];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			for(std::set<atom_index>::const_iterator kt = k_set.begin(); kt != k_set.end(); kt++)
			if(*kt != *it)
			if(*kt != *jt)
			if(topology.path(*jt,*kt))
			for(std::set<atom_index>::const_iterator lt = l_set.begin(); lt != l_set.end(); lt++)
			if(*lt != *it)
			if(*lt != *jt)
			if(*lt != *kt)
			if(topology.path(*kt,*lt))
			{
				potential_list.insert(new harm(p_k,p_phi0,function::ai2ait(*it),function::ai2ait(*jt),function::ai2ait(*kt),function::ai2ait(*lt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}

bool dihedral::harm_fit::read(const std::string & input,const group_t & group,const topology_t& topology,potential_fit_list_t & potential_list)
{
	io::prefix_mark_expression<real_type,false> p_k,p_phi0;
	vector_index_type i_k = -1,i_phi0 = -1;
	std::vector<std::string> ie;
	if(io::match("dihedral harm %1 %2 : %3",input,&p_k,&p_phi0,&ie))
	{
		*p_phi0 = function::deg2rad(*p_phi0);

		if(~p_k) { i_k = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_k); }
		if(~p_phi0) { i_phi0 = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_phi0); }

		atom_index i,j,k,l;
		std::string si,sj,sk,sl;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2,%3 %4",',',*item,&i,&j,&k,&l)) 
		{	
			if((i!=j) && (i!=k) && (i!=l)  &&   (j!=k) && (j!=l)   &&  (k!=l) ) 
				potential_list.insert(new harm_fit(*p_k,i_k,*p_phi0,i_phi0,function::ai2ait(i),function::ai2ait(j),function::ai2ait(k),function::ai2ait(l)));
		} 
		else if (io::match("%1,%2,%3 %4",',',*item,&si,&sj,&sk,&sl))
		{
			const std::set<atom_index> & i_set = group[si];
			const std::set<atom_index> & j_set = group[sj];
			const std::set<atom_index> & k_set = group[sk];
			const std::set<atom_index> & l_set = group[sl];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			for(std::set<atom_index>::const_iterator kt = k_set.begin(); kt != k_set.end(); kt++)
			if(*kt != *it)
			if(*kt != *jt)
			if(topology.path(*jt,*kt))
			for(std::set<atom_index>::const_iterator lt = l_set.begin(); lt != l_set.end(); lt++)
			if(*lt != *it)
			if(*lt != *jt)
			if(*lt != *kt)
			if(topology.path(*kt,*lt))
			{
				potential_list.insert(new harm_fit(*p_k,i_k,*p_phi0,i_phi0,function::ai2ait(*it),function::ai2ait(*jt),function::ai2ait(*kt),function::ai2ait(*lt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}



std::string dihedral::cosf::record_key() const { return "dihedral cos"; }
std::string dihedral::cosf::record_params() const 
{ std::ostringstream os; os << ios_param << p_A << ios_param << function::rad2deg(p_d) << ios_param << p_m; return os.str(); }			
bool dihedral::cosf::read(const std::string & input,const group_t & group,const topology_t& topology,potential_list_t & potential_list)
{
	param_type p_A,p_d,p_m;
	std::vector<std::string> ie;
	if(io::match("dihedral cos %1 %2 %3 : %4",input,&p_A,&p_d,&p_m,&ie))
	{
		p_d = function::deg2rad(p_d);
		atom_index i,j,k,l;
		std::string si,sj,sk,sl;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2,%3,%4",',',*item,&i,&j,&k,&l)) 
		{	
			if((i!=j) && (i!=k) && (i!=l)  &&   (j!=k) && (j!=l)   &&  (k!=l) ) 
			potential_list.insert(new cosf(p_A,p_d,p_m,function::ai2ait(i),function::ai2ait(j),function::ai2ait(k),function::ai2ait(l)));
		} 
		else if (io::match("%1,%2,%3,%4",',',*item,&si,&sj,&sk,&sl))
		{
			const std::set<atom_index> & i_set = group[si];
			const std::set<atom_index> & j_set = group[sj];
			const std::set<atom_index> & k_set = group[sk];
			const std::set<atom_index> & l_set = group[sl];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			for(std::set<atom_index>::const_iterator kt = k_set.begin(); kt != k_set.end(); kt++)
			if(*kt != *it)
			if(*kt != *jt)
			if(topology.path(*jt,*kt))
			for(std::set<atom_index>::const_iterator lt = l_set.begin(); lt != l_set.end(); lt++)
			if(*lt != *it)
			if(*lt != *jt)
			if(*lt != *kt)
			if(topology.path(*kt,*lt))
			{
				potential_list.insert(new cosf(p_A,p_d,p_m,function::ai2ait(*it),function::ai2ait(*jt),function::ai2ait(*kt),function::ai2ait(*lt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}

bool dihedral::cosf_fit::read(const std::string & input,const group_t & group,const topology_t& topology,potential_fit_list_t & potential_list)
{
	io::prefix_mark_expression<real_type,false> p_A,p_d,p_m;
	vector_index_type i_A = -1,i_d = -1 , i_m = -1;
	std::vector<std::string> ie;
	if(io::match("dihedral cos %1 %2 %3 : %4",input,&p_A,&p_d,&p_m,&ie))
	{
		*p_d = function::deg2rad(*p_d);

		if(~p_A) { i_A = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_A); }
		if(~p_d) { i_d = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_d); }
		if(~p_m) { i_m = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_m); }

		atom_index i,j,k,l;
		std::string si,sj,sk,sl;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2,%3 %4",',',*item,&i,&j,&k,&l)) 
		{	
			if((i!=j) && (i!=k) && (i!=l)  &&   (j!=k) && (j!=l)   &&  (k!=l) ) 
				potential_list.insert(new cosf_fit(*p_A,i_A,*p_d,i_d,*p_m,i_m,function::ai2ait(i),function::ai2ait(j),function::ai2ait(k),function::ai2ait(l)));
		} 
		else if (io::match("%1,%2,%3 %4",',',*item,&si,&sj,&sk,&sl))
		{
			const std::set<atom_index> & i_set = group[si];
			const std::set<atom_index> & j_set = group[sj];
			const std::set<atom_index> & k_set = group[sk];
			const std::set<atom_index> & l_set = group[sl];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			for(std::set<atom_index>::const_iterator kt = k_set.begin(); kt != k_set.end(); kt++)
			if(*kt != *it)
			if(*kt != *jt)
			if(topology.path(*jt,*kt))
			for(std::set<atom_index>::const_iterator lt = l_set.begin(); lt != l_set.end(); lt++)
			if(*lt != *it)
			if(*lt != *jt)
			if(*lt != *kt)
			if(topology.path(*kt,*lt))
			{
				potential_list.insert(new cosf_fit(*p_A,i_A,*p_d,i_d,*p_m,i_m,function::ai2ait(*it),function::ai2ait(*jt),function::ai2ait(*kt),function::ai2ait(*lt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}



std::string inversion::potential::record_atom_index() const 
{ std::ostringstream os; os << function::ait2ai(i) << "," << function::ait2ai(j)<< "," << function::ait2ai(k) << "," << function::ait2ai(l); return os.str(); }


std::string inversion::plan::record_key() const { return "inversion plan"; }
std::string inversion::plan::record_params() const 
{ std::ostringstream os; os << ios_param << p_A; return os.str(); }			
bool inversion::plan::read(const std::string & input,const group_t & group,const topology_t& topology,potential_list_t & potential_list)
{
	param_type p_A;
	std::vector<std::string> ie;
	if(io::match("inversion plan %1 : %2",input,&p_A,&ie))
	{
		atom_index i,j,k,l;
		std::string si,sj,sk,sl;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2,%3,%4",',',*item,&i,&j,&k,&l)) 
		{	
			if((i!=j) && (i!=k) && (i!=l)  &&   (j!=k) && (j!=l)   &&  (k!=l) ) 
			potential_list.insert(new plan(p_A,function::ai2ait(i),function::ai2ait(j),function::ai2ait(k),function::ai2ait(l)));
		} 
		else if (io::match("%1,%2,%3,%4",',',*item,&si,&sj,&sk,&sl))
		{
			const std::set<atom_index> & i_set = group[si];
			const std::set<atom_index> & j_set = group[sj];
			const std::set<atom_index> & k_set = group[sk];
			const std::set<atom_index> & l_set = group[sl];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			for(std::set<atom_index>::const_iterator kt = k_set.begin(); kt != k_set.end(); kt++)
			if(*kt != *it)
			if(*kt != *jt)
			if(topology.path(*it,*kt))
			for(std::set<atom_index>::const_iterator lt = l_set.begin(); lt != l_set.end(); lt++)
			if(*lt != *it)
			if(*lt != *jt)
			if(*lt != *kt)
			if(topology.path(*it,*lt))
			{
				potential_list.insert(new plan(p_A,function::ai2ait(*it),function::ai2ait(*jt),function::ai2ait(*kt),function::ai2ait(*lt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}

bool inversion::plan_fit::read(const std::string & input,const group_t & group,const topology_t& topology,potential_fit_list_t & potential_list)
{
	io::prefix_mark_expression<real_type,false> p_A;
	vector_index_type i_A = -1;
	std::vector<std::string> ie;
	if(io::match("inversion plan %1 : %2",input,&p_A,&ie))
	{

		if(~p_A) { i_A = static_cast<vector_index_type>(potential_list.params.size()); potential_list.params.push_back(*p_A); }

		atom_index i,j,k,l;
		std::string si,sj,sk,sl;
		for(std::vector<std::string>::iterator item = ie.begin(); item != ie.end(); item++)
		if(io::match("%1,%2,%3,%4",',',*item,&i,&j,&k,&l)) 
		{	
			if((i!=j) && (i!=k) && (i!=l)  &&   (j!=k) && (j!=l)   &&  (k!=l) ) 
				potential_list.insert(new plan_fit(*p_A,i_A,function::ai2ait(i),function::ai2ait(j),function::ai2ait(k),function::ai2ait(l)));
		} 
		else if (io::match("%1,%2,%3,%4",',',*item,&si,&sj,&sk,&sl))
		{
			const std::set<atom_index> & i_set = group[si];
			const std::set<atom_index> & j_set = group[sj];
			const std::set<atom_index> & k_set = group[sk];
			const std::set<atom_index> & l_set = group[sl];

			for(std::set<atom_index>::const_iterator it = i_set.begin(); it != i_set.end(); it++)
			for(std::set<atom_index>::const_iterator jt = j_set.begin(); jt != j_set.end(); jt++)
			if(*jt != *it)
			if(topology.path(*it,*jt))
			for(std::set<atom_index>::const_iterator kt = k_set.begin(); kt != k_set.end(); kt++)
			if(*kt != *it)
			if(*kt != *jt)
			if(topology.path(*it,*kt))
			for(std::set<atom_index>::const_iterator lt = l_set.begin(); lt != l_set.end(); lt++)
			if(*lt != *it)
			if(*lt != *jt)
			if(*lt != *kt)
			if(topology.path(*it,*lt))
			{
				potential_list.insert(new plan_fit(*p_A,i_A,function::ai2ait(*it),function::ai2ait(*jt),function::ai2ait(*kt),function::ai2ait(*lt)));
			}
		}
		else throw std::runtime_error("Syntax error : "+input);
		return true;
	}
	return false;
}



